/**
 * This version is stamped on May 10, 2016
 *
 * Contact:
 *   Louis-Noel Pouchet <pouchet.ohio-state.edu>
 *   Tomofumi Yuki <tomofumi.yuki.fr>
 *
 * Web address: http://polybench.sourceforge.net
 */
/* deriche.c: this file is part of PolyBench/C */

#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <math.h>

/* Include polybench common header. */
#include <polybench.h>

/* Include benchmark-specific header. */
#include "deriche.h"

/* Array initialization. */
static void init_array(int w, int h, DATA_TYPE *alpha,
                       DATA_TYPE POLYBENCH_2D(imgIn, W, H, w, h),
                       DATA_TYPE POLYBENCH_2D(imgOut, W, H, w, h))
{
    int i, j;

    *alpha = 0.25; // parameter of the filter

    // input should be between 0 and 1 (grayscale image pixel)
    for (i = 0; i < w; i++)
        for (j = 0; j < h; j++)
            imgIn[i][j] = (DATA_TYPE)((313 * i + 991 * j) % 65536) / 65535.0f;
}

/* DCE code. Must scan the entire live-out data.
   Can be used also to check the correctness of the output. */
static void print_array(int w, int h,
                        DATA_TYPE POLYBENCH_2D(imgOut, W, H, w, h))

{
    int i, j;

    POLYBENCH_DUMP_START;
    POLYBENCH_DUMP_BEGIN("imgOut");
    for (i = 0; i < w; i++)
        for (j = 0; j < h; j++)
        {
            if ((i * h + j) % 20 == 0)
                fprintf(POLYBENCH_DUMP_TARGET, "\n");
            fprintf(POLYBENCH_DUMP_TARGET, DATA_PRINTF_MODIFIER, imgOut[i][j]);
        }
    POLYBENCH_DUMP_END("imgOut");
    POLYBENCH_DUMP_FINISH;
}

/* Main computational kernel. The whole function will be timed,
   including the call and return. */
/* Original code provided by Gael Deest */
static void kernel_deriche(int w, int h, DATA_TYPE alpha,
                           DATA_TYPE POLYBENCH_2D(imgIn, W, H, w, h),
                           DATA_TYPE POLYBENCH_2D(imgOut, W, H, w, h),
                           DATA_TYPE POLYBENCH_2D(y1, W, H, w, h),
                           DATA_TYPE POLYBENCH_2D(y2, W, H, w, h))
{
    int i, j;
    DATA_TYPE xm1, tm1, ym1, ym2;
    DATA_TYPE xp1, xp2;
    DATA_TYPE tp1, tp2;
    DATA_TYPE yp1, yp2;

    DATA_TYPE k;
    DATA_TYPE a1, a2, a3, a4, a5, a6, a7, a8;
    DATA_TYPE b1, b2, c1, c2;

    DATA_TYPE two = SCALAR_VAL(2.0);
    DATA_TYPE one = SCALAR_VAL(1.0);
    DATA_TYPE zero = SCALAR_VAL(0.0);
    DATA_TYPE neg_two = SCALAR_VAL(-2.0);

    k = (one - EXP_FUN(-alpha)) * (one - EXP_FUN(-alpha)) / (one + two * alpha * EXP_FUN(-alpha) - EXP_FUN(two * alpha));
    a1 = a5 = k;
    a2 = a6 = k * EXP_FUN(-alpha) * (alpha - one);
    a3 = a7 = k * EXP_FUN(-alpha) * (alpha + one);
    a4 = a8 = -k * EXP_FUN(neg_two * alpha);
    b1 = POW_FUN(two, -alpha);
    b2 = -EXP_FUN(neg_two * alpha);
    c1 = c2 = 1;

#pragma scop
    for (i = 0; i < _PB_W; i++)
    {
        ym1 = zero;
        ym2 = zero;
        xm1 = zero;
        for (j = 0; j < _PB_H; j++)
        {
            y1[i][j] = a1 * imgIn[i][j] + a2 * xm1 + b1 * ym1 + b2 * ym2;
            xm1 = imgIn[i][j];
            ym2 = ym1;
            ym1 = y1[i][j];
        }
    }

    for (i = 0; i < _PB_W; i++)
    {
        yp1 = zero;
        yp2 = zero;
        xp1 = zero;
        xp2 = zero;
        for (j = _PB_H - 1; j >= 0; j--)
        {
            y2[i][j] = a3 * xp1 + a4 * xp2 + b1 * yp1 + b2 * yp2;
            xp2 = xp1;
            xp1 = imgIn[i][j];
            yp2 = yp1;
            yp1 = y2[i][j];
        }
    }

    for (i = 0; i < _PB_W; i++)
        for (j = 0; j < _PB_H; j++)
        {
            imgOut[i][j] = c1 * (y1[i][j] + y2[i][j]);
        }

    for (j = 0; j < _PB_H; j++)
    {
        tm1 = zero;
        ym1 = zero;
        ym2 = zero;
        for (i = 0; i < _PB_W; i++)
        {
            y1[i][j] = a5 * imgOut[i][j] + a6 * tm1 + b1 * ym1 + b2 * ym2;
            tm1 = imgOut[i][j];
            ym2 = ym1;
            ym1 = y1[i][j];
        }
    }

    for (j = 0; j < _PB_H; j++)
    {
        tp1 = zero;
        tp2 = zero;
        yp1 = zero;
        yp2 = zero;
        for (i = _PB_W - 1; i >= 0; i--)
        {
            y2[i][j] = a7 * tp1 + a8 * tp2 + b1 * yp1 + b2 * yp2;
            tp2 = tp1;
            tp1 = imgOut[i][j];
            yp2 = yp1;
            yp1 = y2[i][j];
        }
    }

    for (i = 0; i < _PB_W; i++)
        for (j = 0; j < _PB_H; j++)
        {
            imgOut[i][j] = c2 * (y1[i][j] + y2[i][j]);
#ifndef NO_PENCIL_KILL
            __pencil_kill(y1[i][j], y2[i][j]);
#endif
        }

#pragma endscop
}

int main(int argc, char **argv)
{
    /* Retrieve problem size. */
    int w = W;
    int h = H;

    /* Variable declaration/allocation. */
    DATA_TYPE alpha;
    POLYBENCH_2D_ARRAY_DECL(imgIn, DATA_TYPE, W, H, w, h);
    POLYBENCH_2D_ARRAY_DECL(imgOut, DATA_TYPE, W, H, w, h);
    POLYBENCH_2D_ARRAY_DECL(y1, DATA_TYPE, W, H, w, h);
    POLYBENCH_2D_ARRAY_DECL(y2, DATA_TYPE, W, H, w, h);

    /* Initialize array(s). */
    init_array(w, h, &alpha, POLYBENCH_ARRAY(imgIn), POLYBENCH_ARRAY(imgOut));

    /* Start timer. */
    polybench_start_instruments;

    /* Run kernel. */
    kernel_deriche(w, h, alpha, POLYBENCH_ARRAY(imgIn), POLYBENCH_ARRAY(imgOut), POLYBENCH_ARRAY(y1), POLYBENCH_ARRAY(y2));

    /* Stop and print timer. */
    polybench_stop_instruments;
    polybench_print_instruments;

    /* Prevent dead-code elimination. All live-out data must be printed
       by the function call in argument. */
    polybench_prevent_dce(print_array(w, h, POLYBENCH_ARRAY(imgOut)));

    /* Be clean. */
    POLYBENCH_FREE_ARRAY(imgIn);
    POLYBENCH_FREE_ARRAY(imgOut);
    POLYBENCH_FREE_ARRAY(y1);
    POLYBENCH_FREE_ARRAY(y2);

    return 0;
}
